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ROUTINE  FOR  FINDING  ROOTS  OF  POLYNOMIALS  WITH  REAL  COEFFICIENTS 

ABSTRACT 

A  routine  has  been  developed  which  computes  the  real  and  the  complex  roots 
of  polynomials  with  real  coefficients  up  to  the  tenth  degree.  In  case  of  an 
even  polynomial  it  computes  its  quadratic  factors  and  solves  each  factor  by 
the  quadratic  formula.  In  case  of  an  odd  polynomial  it  computes  one  real  root 
by  locating  it  and  refining  by  Newton' s  algorithm.  Then  it  removes  the 
computed  root  from  the  polynomial  reducing  its  degree  by  one.  The  reduced 
polynomial  is  of  even  degree  and  it  is  solved  by  quadratic  factors.  The 
routine  fails  when  (a)  the  real  roots  are  of  even  multiplicity,  converging  then 
to  wrong  values,  (b)  the  real  roots  are  of  odd  multiplicity,  not  converging  at 
all,  (c)  the  polynomial  is  badly  conditioned  (very  small  changes  in  coefficents 
cause  large  changes  in  the  values  of  the  roots)  converging  then  to  wrong  values. 

The  routine  is  planned  to  compute  frequencies  in  vibrations  problems  which 
involve  complex  roots.  Statistical  evidence  seems  to  indicate  that  only 
polynomials  with  real  roots  can  be  badly  conditioned.  If  this  were  the  case 
the  handicap  (c)  would  be  of  minor  importance  only. 
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A  polynomial  of  n-th  degree  can  be  written  as 

\  n  n-1  n-i 

f(x)  =  x  +  a  x  +  ...  +  a  x  +  an  =  2^  a  x 

i=o 

where  a  =  1,  and  are  real  numbers.  The  roots  of  f(x)  will  be  denoted  by 

*1*  V . v 

The  coefficient  aQ  of  xn  has  been  set  equal  to  one  for  convenience. 

If  n  is  even  all  the  roots  may  be  complex  numbers.  If  n  is  odd  there 
must  be  at  least  one  real  root,  in  which  case  we  start  by  finding  this  real 
root. 

(a)  Locating  One  Real  Root  for  an  Odd  Polynomial. 

We  shall  establish  first  the  bound  for  all  roots,  which  is  given  by  the 
theorem 

UB(upper  bound)  =  Max  ( a.^ )  +1, 

Where  Max  |a^|  is  the  greatest  coefficient  a^,  (i  =  l,2,...,n)  in  absolute 
value .  Thus 

xx  <  UB. 

Let  us  find  two  consecutive  values  of  f(x)  for  read  values  of  x,  say 
f1  =  fCx^)  and  i‘2  =  f(xg),  where  x^  =  -  UB+Axj  Xg  =  -  UB+2Ax,  and  Ax  is  an 
assigned  step  size  of  x,  chosen  sufficiently  small  that  there  is  likely  to 
be  only  one  real  root  between  -  UB  +  kAx  and  -  UB  +  (k  +  l)  Ax  for  any  k.  If 
f^  <  0  and  fg  >  0  then  the  real  root  is  between  x^  and  Xg.  If  f^  and  fg  have 
the  same  sign,  then  compute  f  ,  and  compare  the  signs  of  fg  and  f ^  and  so  on. 
Suppose  that  at  some  value  x^+1,  f  changes  sign,  that  is 

f,  <  0  and  f.  ..  >  0. 
k  k+1 

Then  the  real  root  x^  is  located  between  x^  and 

xk<VVr 
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(b) 


Refining  the  located  Root  by  Newton^  Mat  hod. 

Let  us  call  the  first  approximation  to  the  root  x^ 


x 


(1) 


X 


k+l' 


The  next  approximation  is  then 


.(2)  _  (1) 


=  xv^  -  f (x(l))/f(xv^),  f(xw) 


(D> 


(H^) 


.a) 


and  generally 


x(J+l)  =  x^  -  f(x^)/f'(*^)* 

The  iteration  continues  until  two  consecutive  values  of  becomes  equal 
to  each  other  within  some  small  preassigned  number.  This  establishes  the  first 
real  root  x1  =  . 


(c)  Removing  from  f(x)  the  Known  Real  Root. 

The  next  step  is  to  remove  this  real  root  by  synthetic  decision,  and  find 
the  coefficients  b^  of  the  depressed  polynomial.  f(x)  may  be  written 

foo  =  ai  *n"i  =  (x  -  x )  bi xn-1"1 

i=o  i=o 


Equating  the  coefficients  of  the  same  powers  in  both  members  of  the 
above  identity  gives 

b  =  1. 
o 


b  =a  +  b  .  x, ;  (r  *  l,2,....,n-l). 
r  r  r-1  1  ’  ' 


The  depressed  polynomial  g^x)  = 


n-l-i 


is  of  even  degree. 


(d)  Computing  Quadratic  Factors  of  Even  Polynomials  ■ 

An  even  polynomial  can  alwaya  be  factored  into  real  quadratic  factors. 
Suppose  that 

2 

x  +  p  x  +  q 
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is  a  real  quadratic  factor  of  f(x).  Then  f(x)  may  be  written 

f  (x)  =  22  a  x  =  (x  +  p  x  +  q)  bi  x 

i=o  i=o 

Equating  the  coefficients  of  the  same  powers  in  both  members  of  the 
above  identity  gives  the  coefficients  b1  of  the  reduced  polynomial. 


b  = 
r 


ar  "  pVl 


qb^gj  (r  -  1,2, . , . .  ,n-2) . 


In  the  case  when  the  quadratic 

2 

x  +  p  x  +  q 

is  not  an  exact  factor  of  f(x),  the  factoring  gives 


(x 


+  p  x 


+  q) 


x 


n-2-1 


+  R 


where 

R  =  S1  x  +  Sg. 

Equating  the  coefficients  of  x  gives 


St  =  b  j 
1  n-l 


Pb 


n-l 


Bairstow  devised  an  iteration  scheme  for  improving  the  values  of  p  and 
q,  in  such  a  way  that  the  coefficients  s  and  and  the  remainder  term 
R  approach  zero. 

F^rstow  Algorithm.  Let  p  and  q  be  the  exact  values  of  the  coefficients 
which  will  make  S1  and  Sg  vanish.  The  approximate  coefficients  S1  and  Sg  are 
functions  of  p  and  q  and  a  Taylor  series  expansion  (nonlinear  terms  being 
neglected)  with  p  =  p  +  Ap;  q  =  q  +  Aq  yields 


S1(p,  q)  =  0  =  S1(p,  q)  +  ApCdS-j/dp)  +  AqfdS^/dq) 
Sg(p,  5)  =  0  =  Sg(p,  q)  +  Ap(dSg/dp)  +  Aq(dSg/dq) 


(1) 


Since  S;L  =  bn_1  and  Sg  =  bn  +  pb^ 
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we  have 


dS^/dp  =  dt^/dp;  dS2/dp  =  SbQ/^p  +  p(3bn_1/dp)  +  b^ 

dS-j/dq  =  db^/dq;  dSg/dq  =  3bn/dq  +  pfdb^j/dq) 

Upon  using  the  relation 


b  =  a  -  pb  .-qb  _ 
r  r  r-1  r-2 

and  introducing  the  notation 

-(ab^ap)  =  dr_1 

^r 

we  obtain  the  recursion  formula  for 


dr  =  b^  -  pdr_.^  -  qdj.^;  d_1  =  0;  dQ  =  1;  (r  =  1,  2, . . . .  ,n-2) 
Sb 

A  similar  recursion  formula  for  can  be  obtained,  giving  -(db ^/dq)  =  dy  g 
Hence 

asj^/ap  =  -  dn_2;  asx/aq  =  -  dn_5 

(2) 

ds2/aP  =  -  dn-1  +  pdn_g  +  bn_x;  asg/aq  =  -  dn_g  -  pdn_5 
Substituting  (2)  in  (l)  yields 

(Ap)  dn_2  *  (Aq)  dn_j  -  bn_1 

<a»-l  -  Vl>  +  M  d»-2  ■  b„ 

which,  solved  by  Cramer’s  rule  for  Ap  and  Aq  gives 


Vi 

dn-3 

dn-2 

Vl 

= 

— J 

b 

n 

dn-2 

Aq  = 

<dn-l 

-bn.l>  b„ 

Den 
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where 


Den  = 


n-2 


-7 

n-3 


(d  ,  -b  .)  d  . 
'  n-1  n-1'  n-2 


Since  the  non-linear  terms  in  the  Taylor. expansion  have  been  neglected 
the  new  values  of  p  and  q 


p  =  p  +  Ap;  q  =  q  +  Aq 

will  not  represent  the  true  coefficients  of  the  quadratic  factor.  We  expect 
them  to  be  better  approximations  than  p  and  q  in  the  sense  that  they  will 
make  and  Sg  smaller.  This  procedure  can  be  repeated  until  and  Sg  will 
vanish  within  some  prescribed  small  number. 


Initial  Approximations.  The  convergence  of  the  Bairstow  algorithm  depends 
on  the  initial  approximation  being  close  enough  to  the  true  value.  We  have 
selected  the  following  initial  approximation  to  p  aind  q: 


1  2 

p  =  2(  lan-i l/n)n_1>  1  =  laJ  n 

The  above  approximations  are  averages  derived  from  the  well  known 
properties  of  the  polynomial  coefficients, 


*  n 

an  =  TT 
1=1 


V 


n  2  n 

Vl ' 


This  initial  approximation  gave  convergence  in  very  many  widely  different 
cases  which  were  tested.  The  routine  diverged  or  converged  to  wrong  values 
for  polynomials  with  multiple  real  roots.  Some  of  the  9th  degree  polynomials 
with  all  real  roots  were  badly  conditioned  and  converged  to  wrong  values. 
Polynomials  with  all  or  some  complex  roots  converged  in  all  cases  to  the 
.correct  values. 

It  has  been  found  that  computation  of  quadratic  factors  of  odd 
polynomials  by  the  Bairstow  algorithm,  which  we  use  in  our  routine,  diverged 
in  all  cases  which  were  tested,  even  in  those  where  the  initial  approximations 
were  very  close  to  the  true  values.  This  interesting  fact  has  not  to  my 
knowledge  been  registered  in  the  existing  literature  on  the  subject. 
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3 

High  degree  polynomials  with  coefficients  of  the  order  10  or  higher 
must  have  the  roots  reduced  by  a  constant  factor,  say  10" in  order  to 
sake  the  coefficients  smaller.  Otherwise  the  Intermediate  products  arising 
in  the  computations  may  exceed  in  magnitude  the  capacity  of  a  machine. 

The  appended  flow  chart  and  the  program  in  the  FORAST  language  explain 
the  computational  procedure. 


TADBUSZ  IJ5SER 
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APPENDIX  I 
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FLOW  CHARE 
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FLOW  CHART  (Cont'd) 


Compute  Amax  =  Max  (JaJ,  |a2|  . .  |aJ) 

Compute  Upper  bound  =  UB  =  A  +1. 

max 

Compute  consecutive  values  ffc  =  f  (x^) ;  is  =  1,2,.... 

(initial  value  xq  =  ((K)(Step)  -  UB) . 

Compare  two  consecutive  values  f  ^  and  f. +. 

If  and  have  opposite  sign  then  the  root  x^  is  located,  x^  <  x^  < 

N-l  -  ,  ± 

Find  coefficient  of  the  derivative  polynomial  y  C^x“"  “  . 

(f’(x)  = 

n=o  i=o  * 

ci  =Ai(N-i),(i=1,2...n-1). 

Refine  the  value  of  x^  by  Newton's  method 
Taking  as  initial  approximation  x^  =  x^+1 
x(i+1)  =  x(i)  .  f(x1)/f,(x1)i 

When  X(k+1)  =  T*n*  r  Vi i  n  qmfl  1 1  minjhpr  fhpn 


x  within  small  number  then 


x^  =  x.,;  Print  x^. 


Does  N-l=2 


Find  the  coefficient  of  the  reduced  polynomial  ^  .  T  x 
n-l 

(f(x)  =  (x  -  x2)  T±x  ); 


n-i-1 


+  X;L  (i-l);  Replace  A^^  by  T^. 
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PROGRAM  IN  FORAST  LANGUAGE 


BLOC ( A-A1G )B-B10>D-D10) C-C10 I S-S 1 1 > T-T 10 > *  2 

START  READINFISTEPI*  3 

READ  1 10  ) NOS*  AT  (AH  4 

H®NF/2*  ENTERIWH.FRAIHIWHNIFRANI*  5 

I F(FRAN®0.5) WITHIN! .0000 1 > GOTO (ODDI  %  6 

I F-NOT J  NF-2®0 1  W l TH I N ! . 0000 1 )  GOTOIEVEN1*  7 

ALPl'Al*  BET1®A2*  GOTO! SO) %  g 

ODD  ENTERICVFTOIINFINI*  9 

SET(R«1I*  AM®ABS(A1I*  X0 

1.  I F-ABS I A • I R+ 1 1  *AM)GOTO|2.I*  U 

COUNTIN  I INIRIGOTO! 1. I*  GCTOI3.)  12 

2.  AM®ABS| A.tR+lI I*  13 

COUNT! N) INIRIGOTO! 1. I*  14 

3.  UB®1*AM*  i g 

SUM® 1%  SET ( I ®  X ) *  16 

3.1  SUM®SUM»(-UB)+A,I*  17 

COUNTIN+I I  INI  I IGOTO! 3.1 )%  18 

SUM1®SUM*  SE T ( <« 1 ) *  19 

3.11  ENTERICVITOFIXIKFI*  20 

EX1®KF«STEP-UB  21 

S1®1*  SET! I®ll*  22 

GPS  S » ( I  +  1 1 ®S» I*EXl  +A.I*  23 

COUNT I N  +  l I  IN  I  I  I  GOTO  I  GPS  I  *  24 

SL'M2®St  IN+ll*  25 

I F I  ASS ( SUM 2-SUM 1  I  A® ACS ( SUM2+SUM1 1  I  GOTO! 3.21 »  26 

SUM1"SUM2*  1NT|K°K+1I*  GOTOI3.11I*  27 

3.2  SET! I ® 1 )*  C0®NF*  28 

4.  ENTER(CVITOF) 1  I  IF )*  29 

C. I®A.I»(NF-IFI*  30 

COUNTIN ) INI  I  I  GOTO! 4. )*  31 

SET! J®1 )*  SU1®C0S  32 

GPSD  SO. I J+l l®SD.J»EXl  +C.J*  33 

COUNT  IN ) INI J )  GOTO  I GPSO I  *  34 

SUMO'SD.N*  35 

4.1  EX°EX1  - 1 SUM2 / SUMO  I G  36 

1 F | EX®EX1  I w I THI N I .00001) uOTO I  ROOT  1 1  37 

£X1  “EX*  S1®1*  SET! I ® 1  I  *  SO  1 ®C0*  38 

5.  S.II+ll®S.i®EX+A.I*  SD.|I+1)®SD.I*EX+C.I*  39 

COUNTIN+ll INI  I  I  GOTO  I S. I  *  40 

SUM2“S. (N+l I*  SUMO® SO .NX  GOTOI4.il*  41 

ROOT  1  PR1NT»RR00T*!EXI*  42 

TO® 1*  SET  I I»1 I*  43 

6.  T  » I ®A  » 1 +EX*T  » I  1-1  I  *  44 

COUNT  IN) INI  I ) GOTO  I  6 . I*  45 

MOVE  1 N I  NOS .FROM  I T I  TO  I  A I  *  46 

IF-NOT|NF-l®2 I  WITHIN! .00001  I  GOTO  I  7. I  I  4  7 

ALP1»A1*BET1®A2*NF°NF-1*&OTO(SUIS  48 

7.1  I  NT  I N°N— 1 1  *  NF®NF-1*  GOTUlEVENI  49 

EVEN  ENTERICVFTOIINFINI*  50 

P®ABSIA.(N-1 1  I /NF*Q®ABi I A.N I  *  51 

Y“1/(NF-1I*  LNTERlPO+ERIPIYIALPI*  ALP1®2*ALP*  52 

X«2/NF*  ENTERIPOWCRUIXIBETII*  53 

BE1  INT(N°N  +  1I*  SET ( R®2 1  *  HO'l*  D0®1*  54 

B1®A1-ALP1S  01®B1-ALP1*  55 

BE  B.R®A.R-ALP1»B.(R-1)-BET1»B,(R-2I*  56 

D. R®B.R-ALP1»D, (R-l i-BET1®D.|R-2 I*  57 

COUNT (N)IN(R)GOTCISE)*  58 

1  NT  I  N®N-1 1  *  S 1  ®B  .  I  N-l  I  %  S2®U.N+ALP1®B.(  N-l  I  *  59 

S®A6SIS1I+A3S(S2I*  60 

1 F (G®OI Wl THIN! .0001 1 jOTOISJI *  61 

DCL®0.(N-2l»*2-D.(N-3l®(0.(N-ll-[i,IN-ll  I*  62 

OELAL® I b. (N-l I *0 . ( N-2 l-8.N*0* (N“3 I I /DEL*  63 

DELI3E®lD.lN-2l*B.N-d.(N-l)»(D,(N-ll-B,(N-ll I  I /DEL*  64 

ALP2®ALP1+0ELAL*  BET2®BET1+0ELBE*  65 

ALP1®ALP2*  BET1°BET2*  GCTOIBEII*  66 

SG  DFLTA®ALP1«*2-4»BET1*  67 

IFIDELTA»OlGOTO(8.IS  68 

ALPH®-ALPl/2»  BET®SQRTIDELTA)/2*  69 

X 1®ALPH+BET*  X2®ALPH-BET*  70 

ROOTR  PRINT  *R ROOTS* I X 1  I (X2I*  71 

PRINT*C0EFS*IALP1I I3ET1 I*  72 

I F-NOT INF-2® 01  WITH  IN  I. 0000 II  GOTO  I  9. 1  7  3 

GOTOISTARTI*  74 

3.  ALPH® -ALP  1  /  2S1  BET°SGSTI-0ELTA)/2*  75 

ROOT  I  PRINT*I ROOTS* I ALPHI (BET  I*  76 

PRINT®C0EFS*(ALP1I  (BETH*  77 

I F-NOTI NF-2®0 I W I TH I N( .000C1 1  GOTO  1 9. I  *  78 

GOTOISTARTI*  79 

9.  I F-NOT I NF-2®2 1 Wl THIN! .00001 1  GOTO (9.1 1*  80 

I  NT  I N® N-2 I  *  NF»NF-2*  ALPl'Bl*  BET1«B2*  GOTOIS3I*  81 

9.1  NF®NF-2*  I  NT  I N®N-2 1  *  82 

MOVE  I Nl NOS. F ROM 1 bl I  TUI  A 1  I  *  83 

GOTO! EVEN  I  *  84 

I.NO  GOTOISTARTI*  85 
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